Pseudotime Trajectory Inference with Slingshot
Sedat Kacar
2025-05-05
Welcome to this blog post on performing trajectory inference using the Slingshot package in R. Here, we’ll walk through a full pipeline from preparing your data to visualizing inferred pseudotime and gene dynamics.
1- LOAD LIBRARIES and Prepare Slingshot object
library(Seurat)
library(SingleCellExperiment)
library(ggplot2)
library(dplyr)
SedatFunctions::rds_function(8)
# Assuming you have a Seurat object named 'seurat_obj'
sce_obj <- as.SingleCellExperiment(Endo_sickle)
sce_obj <- sce_obj[rowSums(assay(sce_obj, "counts")) > 0, ]2- CALCULATE TRAJECTORIES using slingshot from the UMAP coordinates.
library(slingshot)
pto <- slingshot(reducedDim(sce_obj,'UMAP'), clusterLabels = sce_obj$seurat_clusters,
start.clus = 3, end.clus = c(4, 5), spar = 1.1)3- Visualizing Slingshot-Inferred Lineage Paths on UMAP Using Minimum Spanning Tree (MST)
slsMST <- slingMST(pto, as.df = TRUE)
scater::plotReducedDim(sce_obj, dimred = 'UMAP', colour_by = 'seurat_clusters') +
ggplot2::geom_point(data = slsMST, aes(x = umap_1, y = umap_2), size = 3, color = "grey15") +
ggplot2::geom_path(data = dplyr::arrange(slsMST, Order),
ggplot2::aes(x = umap_1, y = umap_2, group = Lineage)) +
ggplot2::labs(title = "Slingshot MST on UMAP", x = "UMAP 1", y = "UMAP 2")
4- VISUALIZE ORIGINAL UMAP
hvg_genes <- VariableFeatures(Endo_sickle)
rowData(sce_obj)$highly_variable <- rownames(sce_obj) %in% hvg_genes
sce_obj <- sce_obj[rowSums(assay(sce_obj, "counts")) > 0, ]
scater::plotReducedDim(sce_obj, 'UMAP', colour_by = 'seurat_clusters', text_by = 'seurat_clusters') + theme_bw()
5- Overerlay Inferred Trajectories as Curves
emb <- slingshot::embedCurves(pto, reducedDim(sce_obj,'UMAP'))
emb <- slingshot::slingCurves(emb)
g <- scater::plotReducedDim(sce_obj, 'UMAP', text_by = 'seurat_clusters', colour_by='seurat_clusters')
for (path in emb) {
embedded_curves.df <- data.frame(path$s[path$ord,])
g <- g + ggplot2::geom_path(data=embedded_curves.df, aes(x=umap_1, y=umap_2), linewidth=1.2)
}
print(g)
sce_obj <- scater::runUMAP( sce_obj, dimred = 'PCA',
BPPARAM=BiocParallel::MulticoreParam(6),
ncomponents = 3, name = 'UMAP_3D'
)6- 3D UMAP + Pseudotime Overlay
#3
embedded <- slingshot::embedCurves(sce.sling, "UMAP")
embedded <- slingshot::slingCurves(embedded)
#4
g <- scater::plotReducedDim(sce.sling, 'UMAP', text_by = 'seurat_clusters', colour_by='seurat_clusters')
for (path in embedded) {
embedded_curves.df <- data.frame(path$s[path$ord,])
g <- g + ggplot2::geom_path(data=embedded_curves.df, aes(x=umap_1, y=umap_2), linewidth=1.2)
}
print(g)
#5
sce.sling$Pseudotime <- rowMeans(slingshot::slingPseudotime(sce.sling), na.rm=TRUE)
#6
g <- scater::plotReducedDim(sce.sling, 'UMAP', text_by = 'seurat_clusters', colour_by='Pseudotime')
for (path in embedded) {
embedded_curves.df <- data.frame(path$s[path$ord,])
g <- g + ggplot2::geom_path(data=embedded_curves.df, aes(x=umap_1, y=umap_2), linewidth=1.2)
}
print(g)
#7
embedded <- embedCurves(sce.sling, 'UMAP_3D')7- 3D Visualization of Slingshot-Inferred Lineages Colored by Cell Clusters in UMAP Space
library(dplyr)
library(plotly)
?reducedDim
plot.df <- as.data.frame(reducedDim( sce.sling, 'UMAP_3D' ) )
plot.df$cluster <- sce.sling[['seurat_clusters']]
colnames(plot.df) <- c('UMAP1', 'UMAP2', 'UMAP3', 'cluster')
redDim <- 'UMAP'
p <- plot_ly( ) %>% # colors = trace.colors ) %>%
add_trace( data = plot.df,
x = ~UMAP1, y = ~UMAP2, z = ~UMAP3,
color = ~cluster,
type = 'scatter3d',
mode = 'markers',
marker = list( size = 3, opacity = 1 )
)
for (i in 1:4) {
embedded_curves.tmp <- slingCurves(embedded)[[i]] # only 1 path
embedded_curves.tmp.df <- data.frame(embedded_curves.tmp$s[embedded_curves.tmp$ord,])
p <- p %>%
add_trace( data = embedded_curves.tmp.df,
x = ~UMAP1, y = ~UMAP2, z = ~UMAP3,
type = 'scatter3d', mode = "lines",
name = paste0("Lineage ", i),
line = list( width = 5, opacity = 1 )
)
}
p %>%
layout( scene = list( xaxis = list(title = paste0( redDim, '1') ),
yaxis = list(title = paste0( redDim, '2') ),
zaxis = list(title = paste0( redDim, '3') )
)
)8-Interactive 3D Visualization of Slingshot Trajectories Colored by Pseudotime in UMAP Space
plot.df$Pseudotime <- sce.sling$Pseudotime
p <- plot_ly( ) %>% # colors = trace.colors ) %>%
add_trace( data = plot.df,
x = ~UMAP1, y = ~UMAP2, z = ~UMAP3,
type = 'scatter3d',
mode = 'markers',
showlegend = FALSE,
marker = list( color = ~Pseudotime,
size = 3,
colorscale = 'Viridis', # Viridis color scale
colorbar = list( # Add color scale bar,
title = 'Pseudotime',
x = 1.01, # horizontal position
y = 0.4, # vertical position
len = 0.8, # length (0 to 1, relative to plot height)
thickness = 30 # width in pixels
),
opacity = 1
)
)
for (i in 1:4) {
embedded_curves.tmp <- slingCurves(embedded)[[i]] # only 1 path
embedded_curves.tmp.df <- data.frame(embedded_curves.tmp$s[embedded_curves.tmp$ord,])
p <- p %>%
add_trace( data = embedded_curves.tmp.df,
x = ~UMAP1, y = ~UMAP2, z = ~UMAP3,
type = 'scatter3d', mode = "lines",
name = paste0("Lineage ", i),
line = list( width = 5, opacity = 1 )
)
}
p %>%
layout( scene = list( xaxis = list(title = paste0( redDim, '1') ),
yaxis = list(title = paste0( redDim, '2') ),
zaxis = list(title = paste0( redDim, '3') )
)
)sce.sling <- scater::logNormCounts(sce.sling, assay.type="logcounts")
sce.sling
assayNames(sce.sling)library(TSCAN)
# scran::testPseudotime(...)
pseudo <- TSCAN::testPseudotime( sce.sling,
pseudotime=sce.sling$Pseudotime,
BPPARAM = BiocParallel::MulticoreParam(6)
)9-Filter by Selecting significant genes
library(dplyr)
# ensembldb::filter(., FDR < 0.05)
pseudo.sig <- pseudo %>%
as.data.frame() %>%
dplyr::filter( FDR < 0.05 ) %>%
arrange( FDR ) 10-Visualize top 25 genes on heatmap
topgenes <- rownames(pseudo.sig)[1:50]
scater::plotHeatmap(sce.sling,
order_columns_by='Pseudotime',
colour_columns_by="seurat_clusters",
features = topgenes,
center=TRUE,
scale=TRUE)
11- GENES CHANGING ALONG THE TRAJECTORIES
11a- EXPRESSION OF TEK w/ individual data points
# Example for one gene (e.g., topgenes[1])
gene <- topgenes[39]
df <- data.frame(
pseudotime = sce.sling$Pseudotime,
expression = logcounts(sce.sling)[gene, ],
cluster = sce.sling$seurat_clusters
)
ggplot2::ggplot(df, aes(x = pseudotime, y = expression, color = cluster)) +
ggplot2::geom_point(alpha = 0.5) +
ggplot2::geom_smooth(method = "loess", se = FALSE) +
ggplot2::labs(title = gene, y = "Log-normalized expression")
11b- EXPRESSION OF TEK w/o individual data points
ggplot2::ggplot(df, aes(x = pseudotime, y = expression, color = cluster)) +
ggplot2::geom_smooth(method = "loess", se = FALSE, size = 1.2) +
ggplot2::labs(title = gene, y = "Log-normalized expression")
message("Here is top genes:")## Here is top genes:topgenes## [1] "Nckap5" "Calcrl" "Gm31243" "Bank1"
## [5] "Pax5" "Adgrl3" "Adgre5" "Fendrr"
## [9] "Dennd4a" "Ptprb" "Slfn14" "Tmem100"
## [13] "Slc4a1" "BE692007" "Blnk" "Tspan7"
## [17] "Slc16a10" "Hmbs" "Pim1" "Ldb2"
## [21] "Ikzf3" "Ptprm" "Hmcn1" "Tmem131l"
## [25] "Vwf" "Pde3a" "Epas1" "Btla"
## [29] "5430401H09Rik" "Rgs6" "Kit" "Aff3"
## [33] "Prickle2" "Arl15" "Gpx1" "Blvrb"
## [37] "Prkcb" "Bmpr2" "Ube2o" "Bmp6"
## [41] "Tbcel" "Trim10" "Cdk14" "Tent5c"
## [45] "Eif2ak1" "Ppp1r15a" "Creg1" "Tnfaip2"
## [49] "Hpgd" "Gypa"11c- Plotting Gene Expression of Tek (Tie2) on UMAP
scater::plotReducedDim(sce.sling, dimred = "UMAP", colour_by = topgenes[39])
That’s a wrap! You’ve now built a full Slingshot pipeline, from dimensionality reduction and clustering to trajectory inference, pseudotime visualization, and gene trend analysis.
Happy analyzing!